%%%%%%%%%%%%%%%%%%
% Options return factor model, extended to examine diversification effect
% model: Return = Structural factors (HC, VR, JR)+[HRP, VRP]+ option return [Beta, Idio vol] 
% Liuren Wu, liuren.wu@baruch.cuny.edu
%%%%%%%%%%%%%%%%%%
function [Bi,nobt,Af,Bf,R2s,stats,ccx,ccr,crr,VB] =FunDailyORSestimationDivEffect(t,y,x,HVV,IV,ids,dds,dd,nlag,nx,winth,minN)
tt=dds==dd(t);  sdt=ids(tt);
yt=y(tt);

Nt=sum(tt);
rsv=NaN(Nt,nlag);
jj=1:nlag; nj=length(jj);
Biv=NaN(3,nj);Bi0=ones(3,1)/3;
for j=1:nj
    ll=jj(j);
    ttp=dds==dd(t-ll);sdp=ids(ttp);ytp=y(ttp);
    [~,ind1,ind2]=intersect(sdt,sdp);
    rsv(ind1,j)=ytp(ind2);
    
    HVVp=HVV(ttp,:);
    
    tti=dds==dd(t-ll+1);HVi=HVV(tti,1);sdi=ids(tti);
    [~,ind1,ind2]=intersect(sdi,sdp);
    yi=HVi(ind1);
    Xi=HVVp(ind2,:);
    ii=isfinite(yi+sum(Xi,2)); yi=yi(ii);Xi=Xi(ii,:);
    XX=Xi'*Xi;
    if rank(XX)<3
        Pi=diag(diag(XX))*1e-4;
        Biv(:,j)=(XX+Pi)\(Xi'*yi + Pi*Bi0);
    else
        Biv(:,j)=Xi\yi;
    end
end
Bi=nanmean(Biv,2);

HVVt=HVV(tt,:);FHV=HVVt*Bi;
FE=(IV(tt) -FHV); %future expectation of risk premium

PR=nanmean(rsv(:,2:nlag),2);%moment, past realization

%variance estimates
%rsv Nx12
yr=rsv'; xr=[ones(nlag,1),nanmean(yr,2)]; br=xr\yr;
er=yr-xr*br;
ver=var(er);
nir=find(~isfinite(br(2,:)));
ni=length(nir);
for ii=1:ni
    iii=nir(ii);
    idx=isfinite(yr(:,iii));
    if sum(idx)>5
        br(:,iii)=xr(idx,:)\yr(idx,iii);
        ver(iii)=var(yr(idx,iii)-xr(idx,:)*br(:,iii));
    end
end
svr=br(2,:)';

ivr=ver(:);
crr=corr(nanmean(yr,2), yr,'rows','pairwise');

xt=[svr,ivr, x(tt,:),PR,FE*100]; 

nobt=[Nt,sum(isfinite(yt)),sum(isfinite(xt))];
tvr=nanstd(yr).^2; tvr=tvr(:);

xt=winsorize(xt,winth);
%1. Descriptor stats
Zt=[yt,xt,tvr];
stats=[...
    nanmean(Zt);
    nanstd(Zt,0);
    prctile(Zt,[10,25,50,75,90]');
    skewness(Zt,0,1);
    kurtosis(Zt,0,1)-3;
    100*corr(yt,Zt,'rows','pairwise');
    100*corr(yt,Zt,'rows','pairwise','type','Spearman')];
ccx=corr(Zt,'rows','pairwise');
ccr=corr(Zt,'rows','pairwise','type','Spearman');

%regression
SZt=NaN(Nt,nx); 
indfx=zeros(nx,1);
for j=1:nx
    indf=isfinite(xt(:,j)+yt); Nj=sum(indf);
    if Nj>minN
        indfx(j)=1;
        xtj=xt(indf,j);sxj=std(xtj);mxj=mean(xtj);
       if j<=4
            SZt(indf,j)=xtj./sxj;
        else
            SZt(indf,j)=(xtj-mxj)./sxj;
        end
    end
end

Bf=NaN(nx,1); R2s=NaN(1,1);Af=NaN(1,1);VB=Bf;
indff=find(indfx==1);
if ~isempty(indff)
    Ft = SZt(:,indff);     [Af(1),Bf(indff,1),R2s(1),VB]=FunFMRegression(Ft,yt);
end


end


function [A,B,AR2,vb]=FunFMRegression(Ft,yt)
indf=isfinite(sum(Ft,2)+yt); Nf=sum(indf);

Xtf=[ones(Nf,1),Ft(indf,:)];ytf=yt(indf); 
B=Xtf\ytf; yh=Xtf*B; 
ef=ytf-yh;
nx=size(Xtf,2);
AR2= 1-(sum(ef.^2)./(Nf-nx))/var(ytf);

A=B(1); B=B(2:end);
vb=B.*(cov(Xtf(:,2:end))*B)./var(ytf);


end




